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r~! Context. It is important for the star formation process to understand the collapse of a prestellar dense core. 

Aims. We investigate the effect of the magnetic field during the first collapse up to the formation of the first core, focusing 
i i particularly on the magnetic braking and the launching of outflows. 

Methods. We perform 3D AMR high resolution numerical simulations of a magnetically supercritical collapsing dense core using 
7-H ■ the RAMSES MHD code and develop semi-analytical models that we compare with the numerical results. 

^ ' Results. We study in detail the various profiles within the envelope of the collapsing core for various magnetic field strengths. 
\Q , Even modest values of magnetic field strength modify the collapse significantly. This is largely due to the amplification of the 
radial and toroidal components of the magnetic field by the differential motions within the collapsing core. For a weak magnetic 
intensity corresponding to an initial mass-to-flux over critical mass-to-flux ratio, [i equals to 20, a centrifugally supported 
disk forms. The strong differential rotation triggers the growth of a slowly expanding magnetic tower. For a higher magnetic 
&\ ' field strengths corresponding to /i = 2, the collapse occurs primarily along the field lines, therefore delivering weaker angular 
momentum in the inner part whereas at the same time, strong magnetic braking occurs. As a consequence no centrifugally 
supported disk forms. An outflow is launched from the central thermally supported core. Detailed comparisons with existing 
' analytical predictions indicate that it is magneto-centrifugally driven. 
J> Conclusions. For cores having a mass-to-flux over critical mass-to-flux radio fj, < 5, the magnetic field appears to have a significant 
V ""J ■ impact. The collapsing envelope is denser and flatter than in the hydrodynamical case and no centrifugally supported disk forms. 
' For values fi < 20, the magnetic field drastically modifies the disk evolution. In a companion paper, the influence of the magnetic 
^ , field on the dense core fragmentation is studied. 
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1. Introduction 

Studying the collapse and the fragmentation of a proto- 
stellar molecular dense core is of great relevance for the 
star formation process. While the role of the magnetic field 
has long been suspected (e.g. Shu et al. 1987), it is still a 
disputable issue. 

The first calculations of a collapsing dense core were 
monodimensional and treated ambipolar diffusion (e.g. 
Ciolek & Mouschovias 1994, Basu & Mouschovias 1995). 
Their main goal was to investigate the role of the magnetic 
support in delaying the protostar formation. At about the 
same time, a few attempts were made to calculate the col- 
lapse in 2 or even 3D (Phillips & Monaghan 1985, Fiedler 
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& Mouschovias 1992). In parallel to the numerical efforts, 
various authors have looked for analytical solutions of this 
problem (Galli & Shu 1993ab, Nakamura et al. 1995, Li 
& Shu 1996, Basu 1997, Krasnopolsky & Konigl 2002, 
Hennebelle 2003, Tilley & Pudritz 2003). 

With the increasing computing power and the im- 
provement of the numerical schemes, recent developments 
have been realized and various 2D (Nakamura et al. 1995, 
Tomisaka 1998, Allen et al. 2003) as well as 3D numerical 
calculations have been performed (Hosking & Whitworth 
2004, Machida et al. 2005, 2007, Ziegler 2005, Banerjee & 
Pudritz 2006, Fromang et al. 2006, Price & Bate 2007). 

In these calculations, it has been found that the mag- 
netic field plays a crucial role in the evolution of the col- 
lapsing dense core, in particular in the context of fragmen- 
tation in multiple systems. It has also been found that 
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outflows can be spontaneously launched during the col- 
lapse. These outflows have strong similarities with the one 
studied in many papers either numerically (e.g. Uchida & 
Shibata 1985, Casse & Keppens 2003, Pudritz et al. 2007) 
or analytically (e.g. Blandford & Payne 1982, Pelletier & 
Pudritz 1992, Contopulos & Lovelace 1994, Ferreira 1997). 

Here, we present further 3D numerical calculations 
of a collapsing magnetised dense core. Our main goals 
are to investigate the influence of the magnetic field 
strength on the collapsing envelope, the disk and the out- 
flows. The fragmentation is studied in a companion paper 
(Hennebelle & Tcyssicr 2007, hereafter paper II). In order 
to identify the physical mechanisms at play, we develop 
various analytical approaches that we then compare with 
the numerical solutions. The outline of the paper is as fol- 
low. In the second part, the numerical setup and the initial 
conditions are presented. The third part studies the evo- 
lution of the envelope. For this purpose, semi-analytical 
solutions are obtained and compared with the numeri- 
cal results. The fourth part presents the results for the 
outflows. Comparison with classical analytical results are 
made. In the fifth part, we qualitatively compare our re- 
sults with various observations, focussing particularly on 
the young class source, IRAM04191 (Andre et al. 1999, 
Belloche et al. 2002) The sixth part concludes the paper. 

2. Numerical setup and initial conditions 

We perform 3D numerical simulations using the AMR 
code RAMSES (Teyssier 2002, Fromang et al. 2006). 
RAMSES is based on shock capturing schemes and can 
handle ideal MHD, self-gravity and cooling. It uses the 
constraint transport method to update the magnetic field 
and therefore is preserving the nullity of the divergence 
of the magnetic field. RAMSES has been widely tested 
and gives results comparable to other MHD codes for a 
large set of benchmarks. The AMR scheme offers access 
to the high resolution needed to treat the problem. All the 
calculations performed in the following use the Roe solver. 

The calculations start with initially 64 3 grid cells. As 
the collapse proceeds, new cells are introduced in order 
to ensure that the Jeans length is described everywhere 
with at least 10 cells. Nine levels of AMR are used for a 
total of about 10 6 grid cells and an equivalent numerical 
resolution of 16384 3 . 

Here, we consider simple initial conditions, namely an 
initially uniform sphere in solid body rotation. The mag- 
netic field is initially uniform and parallel to the rotation 
axis. The sphere is embeded into a diffuse medium hun- 
dred times less dense. This makes that the surface of the 
cloud is initially out of pressure equilibrium and there- 
fore expanding. However, since the cloud as a whole, is 
strongly self-gravitating, the collapse is not affected. The 
motivation to start with such simple conditions, some- 
times considered as the standard test case for gravitational 
collapse of dense cores, instead of, for example, with a 
quasi-cquilibrium configuration, is twofold. First, the mag- 
netised collapse has not been widely explored yet and we 



feel it is important at this stage to choose conditions which 
can be easily reproduced by others. Second, unlike in the 
hydrodynamical case, when magnetic field and rotation 
are considered, the age of the structure does influence the 
angular momentum distribution and the structure of the 
field lines. This makes the choice of starting with such a 
structure in near equilibrium also questionable. 

Initially, the ratio of the thermal over gravitational 
energy, a, is about 0.37 whereas the ratio of rotational over 
gravitational energy, /3, is equal to 0.045. These values are 
comparable to standard values quoted in various studies of 
dense cores and are not too far from typical values inferred 
from observations. The cloud temperature is equal to 11 
K. The cloud has a mass of one solar mass, a radius of 
about Ro ~0.016 pc, a density ~ 5 x 10~ 18 g cm -3 giving 
a freefall time, t// ~ 3 x 10 4 years. In the companion paper 
(paper II), anm = 2 perturbation of various amplitudes 
is added. 

The strength of the magnetic field is expressed in 
terms of mass-to-flux over critical mass-to-flux ratio, 
H = (M/$)/(M/$) c , where (M/$) c = ci/3tt x (5/G) 1 / 2 
(Mouschovias & Spitzer 1976) and $ is the magnetic flux, 
ci has been estimated to be about 0.53. The case \i = 1 
corresponds to a cloud just magnetically supported, i.e. 
magnetic forces balance gravitational forces. Various mag- 
netic strengths arc considered in the following, namely, 
fi = 1000 (quasi- hydrodynamical case), \i = 20 (very su- 
percritical cloud), fi = 5 and \i — 2 (highly magnetised 
super critical cloud). 

In order to avoid the formation of a singularity and to 
mimic the fact that at very high density, the dust becomes 
opaque and therefore the gas becomes nearly adiabatic, 
we use a barotropic equation of state: C 2 = (C°) 2 x (1 + 
(p/pc) 4 / 3 ) 1 / 2 , where C s ~ 0.2 km/s is the sound speed 
and p c = 10 g cm~ 3 . Note that Masunaga & Inutsuka 
(2000) demonstrate that this is a good approximation for 
a one solar mass core. 

However, with such an equation of state, the timestep 
in the central part of the cloud becomes so small that it is 
difficult to follow the collapse during a long period of time. 
In order to avoid that problem, we have also performed 
complementary simulations with a critical density p c /l0 = 
10~ 14 g cm" 3 . 

3. Envelope evolution 

In this section, we study the properties of the various fields 
in the collapsing envelope. We first present our notations 
and a simple semi-analytical approach which will be useful 
to understand the simulation results. 

3.1. Analytical model 

Here, we develop a phenomenological model for the pro- 
files of the various fields near the equatorial plane. We 
stress that the main motivation in carrying out such 
analysis is to have models to interpret more accurately 
the complex numerical results. More elaborate models 
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have been developed (e.g. Galli & Shu 1993ab, Li & Shu 
1996, Krasnopolsky & Konigl 2002) assuming mainly self- 
similarity or equilibrium. Since both are restrictive as- 
sumptions and given the complexity of the numerical re- 
sults obtained in the following, it is unclear to which ex- 
tent these models could be used for the purpose of com- 
parison although they undoubtedly provide a sensible hint 
on the physical processes. 

3.1.1. Notation and assumptions 

We consider an initially uniform cloud of mass M , initial 
radius R , in solid body rotation with angular velocity uj 
and threaded by a uniform magnetic field B z parallel to 
the z-axis. 

In the following, we use standard Cartesian coordi- 
nates (x, y, z) and cylindrical coordinates (r, 6, z) therefore 
r = x 1 + y 2 . 

Let h be the scale height of the cloud near the equator, 
we write (see e.g. Li & Shu 1996): 



h(r) = a x r, 
dC 2 



p(r) = 
B z {r) = 



2wGr 2 ' 
H Z C 2 
" VGr ' 



(1) 



where p and B z are the density and z-component of the 
magnetic field near the equator respectively. In the follow- 
ing, it will be assumed that a, d and H z depend weakly 
on r. It is well known that such scaling is a reasonable 
approximation in the envelope during the class-0 phase in 
particular before the rarefaction wave launched at the for- 
mation of the protostar has propagated significantly (Shu 
1977). 

The structures of the radial and azimuthal components 
of the magnetic field are a little more complex. It is well 
known that for symmetry reasons, B r and Be vanish in 
the equatorial plane, z = 0. Their values increase with z 
until they reach their maximum, after which they decrease 
with z. Since here we are interested only in the value near 
the equatorial plane, we write as Krasnopolsky & Konigl 
(2002) 



B r (r, z,t) — H r (r, t) 



C 2 



h(r) 



B e (r,z,t)=H e (r,t) x 



Gr 

c 2 



(2) 



h(r) VGr' 



These two expressions are valid until z ~ h. At higher 
altitude, B r and Be decrease and tend toward their value 
outside the core which in the present simulations is zero. 
Therefore, it is expected that the values of \B r \ and \B$\ 
at a given radius, r, are maximum at the altitude, z ~ 
h(r), max(S r (r, z)) ~ H r C 2 /VGr and max(_B e (r, z)) ~ 
H e C 2 /VGr. 

Thus, in the following, it seems appropriate to dis- 
play the quantities max(B r (r, z))/B z (r, 0) = H r /H z and 
max(B fl (r,*))/fl z (r,0) = H e /H z . 



3.1.2. Axial and radial components of the magnetic 
field 

Since throughout this work, field freezing is assumed, the 
magnetic flux, $, is conserved within the cloud. Therefore: 



/ 



B z x 2irrdr 



B° z x wR 2 



2ttH z (C 2 /VG)R c , (3) 



where R c is the cloud radius at the current time whereas 
i?o is the initial cloud radius. Thus we have: 



H z = (VG/Gf) x B° Z R 2 /(2R C ). 



(4) 



Note that in this expression the cloud radius R c is not 
known. With our choice of initial conditions, R c does not 
evolve much during the class-0 phase and we will assume 
R c ~ i? in the following. This leads to 



Bz(r) = 



B° Z R 



(5) 



The r-component is less straightforward to obtain. Its 
growth is due to the stretching of the field lines by the 
differential motions within the cloud. In the case of a thin 
and isopedic disk, Li & Shu (1997) demonstrated that the 
magnetic flux and gravitational potential are proportional 
through the cloud allowing one to compute all components 
of the magnetic field once the gravitational potential is 
known. Krasnopolsky & Konigl (2002) have assumed that 
B r is simply proportional to the magnetic flux. Since the 
B r component appears difficult to predict quantitatively, 
we simply write 

H r = r,H z (6) 
and the value of 77 will be estimated from the simulation. 

3.1.3. Density field 

In order to estimate the density, we write axial and ra- 
dial equilibrium conditions. Although the cloud is not ex- 
actly in equilibrium since it is collapsing, such assumptions 
lead nevertheless to reasonable estimates of the density as 
long as the collapse is not strongly triggered (Shu 1977, 
Hennebelle et al. 2003). 

The equilibrium along the z-axis, neglecting the az- 
imuthal component of the magnetic field and the tension 
term B z d r B r , is: 



f;VM (>0-.r- 0, [ ^ 



0. 



(7) 



where ip is the gravitational potential. Integrating once, 
we obtain (using d z ip ~ —AnGp): 



(8) 



where K(r) is a function of r. Evaluating K at z — and 
at z = h, and using the expressions stated by Eqs. (1) 
and (4), we get 
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d~d 2 a 2 + ^-H 2 , 

4 z ' 



(9) 



where we have also used the approximation d z ip ~ 
-AirGph. 

The equilibrium along the radial direction is (neglect- 
ing again the influence of Bg) 

-C 2 s d rP + ^-(-B z d r B z + B z 8 z B r ) + pd r ^ ~ 0. (10) 

47T 



Thus we obtain, with Eqs. (1) and (2) 
H 2 

d + -^(1 + rj/a) ~ ad 2 , 



(11) 



where the gravitational force d r ip has been assumed to be 
d r tp ~ -GM/r 2 with M ~ J 2nr x 2h(r)pdr. 

H z being known from Eq. (4), Eqs. (9) and (11) can 
be solved numerically once r] is estimated from the simula- 
tion, to provide the values of d and a. For the case H z = 0, 
we have a = d = 1, i.e. the structure of the cloud is not 
modified by the magnetic field and therefore the density 
is the Singular Isothermal Sphere (SIS) density (since the 
analytical model does not consider the effect of rotation). 



Table 1. Values for the parameters, r\ is estimated from the 
numerical simulation. H z , a and d are obtained from the ana- 
lytical model. 



A* 


H z 


V 


a 


d 


20 


0.41 


5 


0.48 


2.48 


5 


1.64 


3 


0.2 


10.29 


2 


4.10 


2 


0.12 


29.37 



Gathering Eqs. (1), (2) and (12), we get 
H z H e (t) 



d 7 L 



K x 



ad 



(16) 



where K = C 2 R a /(2GM ). 

Finally, the induction equation together with Eqs. (1) 
and (2) leads to 



d T H e 



-V 



L(t)H z 
f(t) 2 



(17) 



Once a and d are known, Eqs. (14), (16) and (17) can 
be integrated with time to obtain the particle momentum. 
In the following, we use these equations to directly com- 
pare with the numerical results. 



3.1.4. Azimuthal magnetic field and rotation 

The azimuthal component of the magnetic field, as well 
as the rotation are more difficult to obtain. In order to 
do so, we adopt a Lagrangian approach, i.e. we follow the 
fluid particle and compute its momentum and azimuthal 
magnetic field along time. For this purpose, we simply use 
the fluid equations with density and poloidal field given 
as described above. To use dimensionless quantities, we 
define 

r = r/R Q ,M = M(r)/M Q ,t = tx ^GMq/RI (12) 

To compute the position of the fluid particle, we simply 
write (neglecting the thermal pressure) 

(13) 

with V r = dr/dt. In this expression, M(r ) is the mass 
of the cloud within a radius r and G c fj is the effective 
gravitational constant G c e = G x (1 — It will be 

assumed that M remains constant during the collapse, i.e. 
we do not consider any accretion which may arise along 
the pole. Thus, we obtain 



d T V r 



r(t) 2 r(t) 3 



(14) 



where L = rVg is the momentum of the fluid particle. 
The momentum equation is 

MrVe) = ^(B r d r (rB e ) + rB z d z Bg) (15) 



3.2. Cloud radial profiles 

Figures 1-4 show the density, radial velocity, rotation ve- 
locity and z-component of the magnetic field in the equato- 
rial plane (variations along the z-axis are shown in Sect. 4 ) 
as a function of radius for various magnetic field strengths. 
They also display the largest value of radial and azimuthal 
components of the magnetic field at a given radius. These 
values are obtained by taking the largest values along the 
z-axis at each radius. Note that as recalled in the previous 
section, B r and Bg vanish in the equatorial plane z = 0. 
Therefore, the maximum value of B r (r, z)/B z (r, 0) at a 
given r is plotted. Four snapshots are displayed. The first 
one is representative of the prestellar phase and is about 
0.06-0.08xTff before density reaches the critical density, 
p c , the second one is near the time at which the density 
reaches p c whereas the third and fourth ones show latter 
evolution. The two straight solid lines in the density plots 
show the density of the singular isothermal sphere (lower 
lines) and the density of the analytical model stated by 
Eq. (1) (upper lines). Note that in the hydrodynamical 
case, the two straight lines are indistinguishable. Table 1 
gives the values of the parameters, (i, H z , rj, a and d. The 
straight solid lines in the B z plots show the analytical 
estimate of B z stated by Eq. (1) and (4). 

3.2.1. Quasi-hydrodynamical case 

Figure 1 shows results for p, = 1000, i.e. quasi- 
hydrodynamical case. The density is slightly stiffer than 
r~ 2 in the outer part where it is a little denser than the 
SIS. This is due to the rotation and to the fact that the 
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H= 1000 , /?= 0.045 




-4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 
log(r) (pc) 



Fig. 1. Case fi = 1000. Density (left-top), radial (left-middle) 
and azimuthal velocity (left-bottom) and z-component of mag- 
netic field (right-top) in the equatorial plane at different 
times. Largest values of the radial (right-middle) and az- 
imuthal (right-bottom) magnetic components at a given ra- 
dius are also given. For convenience, max(B r (r, z))/B z (r, 0) 
and max(£>e(r, z))/B z (r, 0) are given as a function of r. The 
various straight lines show analytical values (see text). 



cloud is collapsing and therefore not in equilibrium (Shu 
1977, Hennebelle et al. 2003, 2004a). In the inner part of 
the envelope the ratio of density over SIS density increases 
even more with radius. This is due to the rotation velocity 



fi= 20 , (3= 0.045 




-4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 
log(r) (pc) 

Fig. 2. Same as Fig. 1 for case /i = 20. 

which increases with r (Ulrich 1976). Note that a better 
agreement between analytical and numerical estimate can 
be obtained by taking into account the influence of rota- 
tion in the former (see e.g. Hennebelle et al. 2004a). Two 
accretion shocks are visible in the radial velocity plot. The 
first one which is located at r = 10~ 3 pc shows the edge of 
the centrifugally supported disk. The second one, located 
at r = 10~ 4 pc, shows the edge of the thermally supported 
core. Although for this case, the magnetic field has almost 
no influence on the gas dynamics, it is worth studying the 
spatial dependence of the three components. The B z com- 
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fi= 5 , /?= 0.045 




> -1.2 r V 




-4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 
log(r) (pc) 

Fig. 3. Same as Fig. 1 for case fj, = 5. 

ponent in the envelope appears to be reasonably close to 
the analytical estimate stated by Eq. (1), the discrepancy 
being due to the fact that p is stiffer than r~ 2 because 
of rotation. The B r and Bg components which vanish ini- 
tially, have slightly different behaviour. They grow with 
time and reach values of the order of B z in the outer part 
of the envelope down to values roughly 10-100 times larger 
than B z at the edge of the disk. Inside the centrifugally 
supported disk these values increase further up to values 
as high as about 10 3 . It should be noted however that here 
we are plotting the maximum values of B r and Bg at a 



^=2,0= 0.045 




-4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 
log(r) (pc) 



Fig. 4. Same as Fig. 1 for case /i = 2. 



given radius. Since in the case jj, = 1000, the disk quickly 
fragments (see paper II), B r and Bg fluctuate significantly 
and therefore the high values reached in the inner parts are 
much higher than the mean values of B r and Bg (sec paper 
II for an estimate of their mean values in the disk). Note 
also that the increase of max(B r )/ B z and max(Bg)/ B z at 
large radius (r > 10 -2 pc) is simply due to the decrease 
of B z in the external medium surrounding the cloud. 
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3.2.2. Weak field case 

Figure 2 shows results for /i = 20, i.e. weakly magnetised 
case. As expected, since the magnetic field is weak, the 
density, radial and azimuthal velocity fields are very close 
to those obtained in the previous case. B z is much larger 
than in the case /Lt = 1000. As for the previous case, the 
values of max (£?,.)/ B z and maji(Bg)/ B z increase gradu- 
ally in the envelope. They reach values of roughly 10 at 
the edge of the disk. This indicates that the differential 
motions within the cloud are less important in this case 
because of the influence of the Lorentz force. As in the hy- 
drodynamical case, a centrifugally supported disk formed 
as well as two accretion shocks. Note again that the large 
fluctuations of B r and Bg within the disk are due to the 
display of the maximum of B r and Be at a given radius. As 
shown in paper II, some symmetry breaking is occurring 
in the disk which generates strong fluctuations. 
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3.2.3. Intermediate field cases 

Figures 3 and 4 respectively show results for fi = 5 and 
fi = 2, i.e. intermediate and strongly magnetised super- 
critical cases. The density and velocity fields are signifi- 
cantly different from the two preceding cases. The equa- 
torial density is roughly 10 to 30 times the density of the 
SIS and is in good agreement with the analytical estimate 
stated by Eq. (1). This is mainly due to the magnetic pres- 
sure (due to B r ) which compresses the gas along the z-axis. 
In the outer part, the radial velocities are smaller than in 
the weak field cases. This is due to the influence of the 
Lorentz force which supports the cloud. On the contrary, 
in the inner part, the radial velocities are larger than in 
the weak field cases. This is because, since the rotation is 
much smaller than in the weak field case, the centrifugal 
support is much weaker. In the case ji = 5, a weak local 
maximum, due to the centrifugal force is nevertheless still 
present at r ~ 2 — 3 x 10~ 4 pc. However, unlike in the 
cases [i = 1000 and 20, the radial velocity does not van- 
ish except in the center. This indicates that there is no 
real centrifugally supported disk in this case. For /U = 2, 
only the shock on the thermally supported core remains, 
indicating that the centrifugal force is not able to stop the 
gas. The reason for lower angular momentum will be ana- 
lyzed in the next section. We note that similar conclusions 
have been recently reached by Mellon & Li (2007). It is 
also apparent that the shape of the rotation velocity is 
flatter when /j, is smaller: the rotation curve stays roughly 
constant until much smaller radii. 

The z-component of the magnetic field is very close to 
the analytical estimate in the envelope of the core. The 
value of max(B r )/ B z is about 2 at the edge of the core 
for /z = 5 and about 1 for fj, = 2. It gradually increases 
inwards and reaches values about 2-3 times larger in the 
inner part. The values of ma,x(Bg)/ B z are typically 1.5 to 
2 times smaller than max(5 r )/ B z for /j, = 5 and about 3 
times for fi = 2. 



Fig. 5. Fraction of cloud mass, f m , enclosed within a cylinder 
of radius, r/, as a function of 77. Upper panel is case p = 20. 
Lower panel is case /i = 2. 

Altogether these results illustrate that even for low to 
intermediate values of the magnetic strength, magnetic 
field can have a drastic influence on the cloud evolution 
as well as on the disk formation. This is due to the fact 
that the radial and toroidal components of the magnetic 
field, which vanish initially, are strongly amplified during 
the collapse by the differential motions. This makes that 
the radial component B r docs not increase linearly with 
the initial magnetic field strength since the field is easier 
to stretch when it is initially weaker. 

It is worth stressing that such values of \x in the range 
5 — 2 are compatible with the more pessimistic estimates 
derived from measurements of the magnetic intensity in 
the dense cores (Crutcher 1999, Crutcher & Troland 2000, 
Crutcher et al. 2004). Since we find that dense cores having 
/U smaller than ~5 are qualitatively different from the hy- 
drodynamical cores, we conclude that magnetic fields are 
playing an important role in the collapse of dense cores 
and therefore in the star formation process. 

3.3. Angular momentum evolution 

Here, we further study the radial distribution of angular 
momentum. In particular, we investigate the physical ori- 
gin of smaller rotation velocities in the intermediate and 
strong field cases. 

3.3.1. Mass and angular momentum distribution 

For this purpose, we plot the fraction of mass, f m , enclosed 
inside cylinders of various radii for the cases /Lt = 20 and 
/i = 2 in Fig. 5. Note that the first two times correspond 
to a critical density equal to p c , whereas for the three oth- 
ers, the critical density is p c /10. As can be seen, a good 
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Fig. 6. Mean specific angular momentum enclosed within a 
cylinder of radius, r/, containing a cloud mass fraction f m , as 
a function of f m . Upper panel is case fi — 20. Lower panel is 
case [i — 2. 
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Fig. 7. Case /j, — 20. Specific angular momentum of the fluid 
particle located at cloud radius, 77, where r/ is the radius of 
the cylinder enclosing a constant mass fraction / m . Each point 
corresponding to a given time, this diagram shows the time 
evolution of the angular momentum of the fluid particle as the 
collapse proceeds. The solid lines correspond to the analytical 
model presented in Sect. 3 (see text). The dashed lines show 
the radius of the magnetized area which surrounds the ther- 
mally supported core, below which the analytical model is not 
appropriate. 



agreement is obtained between the second and the third 
times which are close in time, showing that varying the 
critical density does not affect significantly the envelope 
evolution. We define the radius r/ of the cylinder contain- 
ing a constant mass fraction f m . rf decreases with time 
as the collapse proceeds. 

Comparison between the 2 panels of Fig. 5 reveals that 
the mass distribution as a function of radius is significantly 
different in the two cases /x = 20 and pL = 2. In particu- 
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Fig. 8. Same as Fig. 7 for case /i = 2. 



lar, the mass fraction enclosed within a cylinder of radius 
~3x 10~ 3 pc is roughly fifty percent higher for [i = 20 
than for [i = 2. This indicates that the collapse arises in 
different ways for these two cases. In fact, the collapse is 
more spherical for /j, = 20 than for /j, = 2. In the latter, 
since the field is strong, the collapse first proceeds along 
the field lines implying that the material which consti- 
tutes the central core and the disk was originally located 
closer to the rotation axis than for the case \i = 20. Since 
material close to the rotation axis has a lower angular mo- 
mentum than the gas located further away, we believe that 
this is one of the reason of lower angular momentum in 
the cloud central parts in the case fi = 2 than in the case 
/i = 20. 

Figure 6 displays the distribution of mean specific an- 
gular momentum within a cylinder enclosing a mass frac- 
tion f m as a function of f m . It shows that the specific 
angular momentum for the first two times are both pro- 
portional to the mass fraction and that the cloud with 
fj, = 20 has a specific angular momentum which is only 
20% higher than for the case \x = 2. This difference which 
is attributable to the magnetic braking, shows that mag- 
netic braking plays only a minor role during the early 
parts of the collapse. Fig. 6 together with Fig. 5 also 
demonstrates that the total angular momentum in the 
case 11 = 20 is higher in the internal part of the cloud 
than in the case fi = 2. The subsequent times shown in 
these figures reveal that the specific angular momentum 
stays roughly constant indicating that magnetic braking 
remains rather inefficient in the case n = 20. On the con- 
trary, the subsequent times displayed in the second panel 
of Fig. 5 show that the specific angular momentum de- 
creases drastically in the inner part of the cloud. This is 
most likely due to the magnetic braking. We stress how- 
ever, that at time t = 1.53r//, the angular momentum 
has decreased significantly only in the very inner part cor- 
responding to f m < 0.15. With Fig. 6, we see that this 
corresponds to radii smaller than 10 -3 pc. Therefore, the 
small rotation velocities seen in Fig. 4 are largely due to 
the collapse proceeding first along the field lines. At this 
time, magnetic braking has efficiently reduced the gas an- 
gular momentum in the very inner parts only. 
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3.3.2. Comparison between analytical and numerical 
results 

In order to compare the numerical results with the ana- 
lytical model and to confirm that the decrease of angular 
momentum seen for fj, = 2 is due to magnetic braking, 
we have calculated the specific angular momentum of the 
fluid particle located at the radius, r/, where r/ is the ra- 
dius of the cylinder which contains a cloud mass fraction 
f m - Indeed, Fig. 6 shows that the mean angular momen- 
tum enclosed within a cylinder of radius Tf does not vary 
much with time (except for the 2 last times displayed for 
/i = 2). This implies that the mass enclosed within radius 
rf, does not vary significantly along time. Therefore, the 
selected fluid particle located at r/ should remain nearly 
the same. Consequently, any angular momentum variation 
is attributable to magnetic braking. Figures 7 and 8 show, 
for different values of f m , the specific angular momentum 
of the fluid particle as a function of radius at ten different 
times. It also displays analytical curves performed with 
the model presented in Sect. 3.1.4. To obtain these curves, 
we start with values of r(t) and L(t) corresponding to the 
first point shown in each panel of Figs. 7 and 8, and we 
integrate Eqs. (14), (16) and (17) using the values of a 
and d provided by Eqs. (9) and (11). Note that in order 
to mimic the growth of B r and the fact that it is initially 
zero, we use r] = r] x (1 — r(t)/r(0j) in Eq. (17). Since as 
shown in Fig. 2, the value of 770 varies through the cloud, 
we run three models for 770=1, 1.5 and 2. The top curves 
of Fig. 8 correspond to t]o = 1 whereas the bottom curves 
correspond to r] = 2. 

The ten times represented in the case \x — 20 cor- 
respond to 1, 1.15, 1.19, 1.2, 1.26, 1.35, 1.46, 1.54, 1.6, 
1.63t//, whereas for \i = 2, they correspond to 1.21, 1.41, 
1.53, 1.52, 1.55, 1.57, 1.60, 1.62, 1.64, 1.68t//. Note that 
for both cases, the first three times have been obtained 
with the standard critical density whereas the seven oth- 
ers have been obtained with the critical density p c /10. The 
good continuity shows that the results are not affected by 
the thermally supported core (except maybe for the last 
times). 

In the case \x = 20, there is, as expected, hardly any 
variation of angular momentum. The only variation occurs 
for f m = 0.1 at radius smaller than r = 3 x 10~ 4 pc, i.e. af- 
ter the fluid particle has reached the central thermally sup- 
ported core. In the case fj, = 2, magnetic braking is much 
more effective. A significant loss of angular momentum is 
observed for f m = 0.1, f m = 0.2 and f m = 0.3. In each 
case, the analytical fit is in reasonable agreement with the 
numerical value until the fluid particle reaches a strongly 
magnetised area surrounding the thermally supported core 
where the analytical solution becomes inappropriate. This 
agreement shows that the analytical model is reasonably 
accurate and that magnetic braking is responsible for the 
angular momentum decrease. Depending on the value of 
f m , the angular momentum decrease during the collapsing 
phase can be larger, comparable or smaller than the angu- 
lar momentum decrease once the fluid particle has reached 



the magnetised area which surrounds the thermally sup- 
ported core. Note that the size of this area increases with 
time due to accumulation of magnetic flux. This is why 
fluid particles corresponding to bigger f m reach it at larger 
radii. 

To summarize, we can say that for low magnetic 
strengths, magnetic braking is too small to play a signifi- 
cant role in the envelope. In the case of strong fields, the 
collapse first occurs along the field lines therefore deliver- 
ing a low angular momentum in the inner region. At the 
same time, magnetic braking reduces the angular momen- 
tum of the collapsing envelope. Finally, strong magnetic 
braking occurs in the surrounding of the thermally sup- 
ported core which is highly magnetised. 

4. Outflows 

It is now well known that accretion is often, maybe always, 
associated with ejection processes. In the context of star 
formation, molecular outflows as well as jets have been 
extensively observed (see e.g. Bally et al. 2007 for a recent 
review). While jets may have velocities as large as several 
hundred km s _1 , the bulk of the millimeter wavelength CO 
emission, tends to have velocities of only a few to ten km 
s _1 . In the following, we call outflows, outward motions 
with velocities larger than 1 km s _1 . 

In this section, we study the various outflow motions 
obtained in these simulations. We first give a basic de- 
scription of the weak and strong field cases. Then, a more 
detailed analysis is presented for each of these two cases. 
Finally, we show for both cases mass and angular momen- 
tum fluxes for long time evolution. 

4.1. Weak field 

Here we describe results for the case [i = 20. 
4.1.1. Basic description 

Figure 9 shows the density and velocity fields in the xz 
plane for /i = 20 at three times. A complex expanding 
structure forms around the center. As will be seen later, 
it is somehow similar to the magnetic tower investigated 
by Lyndcn-Bcll (1996) and in the following we use this 
terminology. The first snapshot shows that this structure 
encompasses the centrifugally supported disk. As a conse- 
quence the accretion shock, which occurs near the equa- 
torial plane in the hydrodynamical case, is located further 
away at the edge of the tower. At this stage, the tower is 
uniformly slowly expanding (see next section for quantita- 
tive estimates). The second snapshot shows that a faster 
outflow appears along the pole. It is clearly starting from 
the central thermally supported core. The velocity is all 
the way almost parallel to the z-axis. Since this outflow is 
faster than the surrounding slowly expanding tower, the 
shape of the structure becomes gradually more complex 
and mainly composed of two distinct regions, the faster 
flow and a slower magnetic tower. The third snapshot 
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Fig. 9. Case pi = 20. Density and velocity fields in the xz plane. 



shows that this structure is maintained at later times with- 
out much change for the central flow whereas the tower 
keeps expanding. At the edges of the structure, near the 
equatorial plane, slow recirculation flows develop. 

It should be noted at this stage that the thermal struc- 
ture of the protostar is not correctly treated in this pa- 
per. In particular the second collapse is not considered 
here (Masunaga & Inutsuka 2000, Machida et al. 2007). 
Thus the central outflow may have a different structure 
in a more realistic simulation. Indeed, Banerjee & Pudritz 
(2006) and Machida et al. (2006, 2007) found that a fast 
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Fig. 10. Structure of the azimuthally averaged magnetic field 
in the model having /i = 20 at time t — 1.2287r//. The solid 
lines displays the poloidal magnetic field lines. They are over- 
plotted on a snapshot of the toroidal magnetic field strength. 



outflow, maybe a jet, having velocities around 30 km/s 
develops during the formation of the protostar. 

Figure 10 shows the structure of the magnetic field at 
time t — 1.2287r//. The magnetic field is decomposed into 
its toroidal and poloidal parts, Bg and B p = (B r ,0, B z ). 
The strength of the former is shown using the colorscale 
snapshots while the poloidal magnetic field lines are rep- 
resented using the solid white lines. The structure of the 
magnetic field appears to be complex. The field lines are 
strongly bent and twisted in the inner central regions 
(x/R c < 0.1, \z/R c \ < 0.1) whereas they are almost 
straight in the outer part. In the same way, the toroidal 
component is 2 to 3 order of magnitude higher in the in- 
ner part than in the external part. This strongly suggests 
that the growth of the tower as well as the outflow are as- 
sociated with the strong wrapping of the field lines. This 
effect is quantified in the following section. 



4.1.2. Quantitative estimates 

Figure 11 shows the density, axial velocity, rotation veloc- 
ity and toroidal component of the magnetic field along 
the z-axis for four times at x/Rq = 0.02 and y = Q. 
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Fig. 11. Cut along the z-axis at x/R c — 0.02 and y — at 
four different times. Density, vertical component of the velocity, 
rotation velocity and toroidal component of the magnetic field 
are shown. 



The first and second times (respectively solid and dot- 
ted lines) show that the central density is increasing due 
to the rapid accretion. Similarly, the angular momentum 
increases. The toroidal component of the magnetic field 
grows rapidly and is about 5 times larger at the second 
time than at the first. This induces a slow expansion at 
about 0.3-0.5 km/s. 

The third time shows that the tower keeps expanding 
with about the same velocity and that the toroidal com- 
ponent of the magnetic field does not grow in intensity 
and saturates, forming a plateau (except close to z = 0) 
with slowing decreasing value at high z. The total toroidal 



magnetic flux inside the structure increases since the size 
of this plateau increases. 

To characterize the dynamical state of the tower, we 
estimate the thermal and magnetic pressure as well as the 
gravitational potential at z ~ 3 x 10~ 3 pc, i.e. close to 
the edge of the tower at time, t — 1.15r//. The density is 
about 10~ 15 g cm~ 3 , giving a thermal pressure of about 
5 x 10 -7 erg cm~ 3 . The toroidal component of the mag- 
netic field is about 10 4 p,G giving a magnetic pressure of 
about Bg /8tt ~ 4 x 10~ 6 erg cm~ 3 . The gravitational force 
is less straightforward to estimate. By the time we are con- 
sidering, the mass denser than ~ 10~ 15 g cm~ 3 is of the 
order of ~ 0.1M S . Thus, the potential energy is of the 
order of ~ pGM/z ~ 10 -5 erg cm -3 . Therefore, we con- 
clude that by the time t — 1.15t//, the magnetic tower is 
largely dominated by the magnetic and the gravitational 
energies. At later times, as the expansion proceeds, the 
gravitational energy will eventually become negligible. 

To assess that the expansion of the tower is indeed due 
to the growth of the toroidal magnetic field, we consider 
pressure equilibrium at the edge of the tower where we 
have 



B 2 /8n = p int V* { , 



(18) 



since the external pressure is dominated by the ram pres- 
sure, pinfV^f, exerted at the accretion shock. 

The flux per unit length of the toroidal magnetic field, 
$0, is given by 



$ 



Bgdz ~ Be x I, 



(19) 



I being the height of the magnetic tower. 

Integrating the induction equation along z, &g is also 
given by 



$ e ~ B z Vg x t, 



(20) 



where B z and Vg are to be taken in the equatorial plane. 
Therefore, we obtain: 



Vtowcr — 



B, 



inf 



(21) 



This expression is somehow similar to some of the ex- 
pressions obtained by Lyndcn-Bcll (1996, 2003) although 
his analysis is more sophisticated since the explicit value 
of the tower radius is taken into account. With B z ~ 
3 x 10 3 pG (obtained from Fig. 2), Vg ~ 0.9 km/s, V int ~ 
1.3 km/s and p m i — 3 x 10~ 17 g cm~ 3 (obtained from 
Fig. 11 either at z = or at z = 9 x 10 -4 pc), we ob- 
tain: Vt owor ~ 0.58 km/s. This value is comparable with 
the value of V z ~ 0.45 km/s at time t = 1.15r// and 
z = 8x 10~ 4 pc within about 25%. The difference is prob- 
ably due to the assumption of constant Bg in the tower. 
Note that this simple estimate does not take into account 
gravity. In order to investigate its influence, an analytical 
model for the expansion of the magnetic tower, is devel- 
oped in the Appendix. Indeed, the model shows that the 
growth of the transverse magnetic field which is induced 
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by the gradient of transverse velocity along the z-axis, 
triggers the expansion of a self-gravitating layer in a very 
similar way to what is observed in the simulation. 

The last time in Fig. 11 shows that the z- velocity in- 
creases significantly and reaches values of about 1.2 km/s. 
This is due to the central outflow which presents higher 
velocities. At this stage the velocities of the tower and the 
flow are difficult to distinguish. The fourth time also re- 
veals that angular momentum as well as mass have been 
removed probably by the outflow between z — 3 x 10 
pc. 



and z = 3 x 1CT 3 



4.2. Strong field: magneto-centrifugal ejection 

Here we describe results for the case /i = 2. 

4.2.1. Basic features 

In the case \i = 2, a collimated outflow developed quickly, 
as seen in Fig. 12. The first time displays the density and 
velocity fields just before the outflow is launched. The sec- 
ond time shows the early phase of the flow whereas the 
third time shows a more advanced phase after which the 
flow characteristics do not evolve much (see next section) . 

The morphology of the flow is quite different from what 
is obtained in the previous case. In particular, there is no 
slow magnetic tower as in the previous case. This is be- 
cause, as discussed in Sect. 3, there is no centrifugally 
supported disk, instead a collapsing magnetic pseudo-disk 
forms. The outflows seem to emerge from the central ther- 
mally supported core with an angle of about 40 — 45 de- 
grees with respect to the z-axis and quickly recollimates. 
Figure 13 shows the structure of the magnetic field lines 
and strength of the toroidal magnetic field. The poloidal 
magnetic field is seen to be mostly vertical, particularly 
away from the equatorial plane. Close to the equatorial 
plane, the magnetic field lines are significantly inclined 
because of the inflowing fluid motions. 

This is qualitatively in good agreement with the now 
classical model of the magneto-centrifugal ejection first 
described by Blandford & Payne (1982) and obtained in 
many simulations of magnetised disks (e.g. Pudritz et al. 
2006). In the following section, a more quantitative anal- 
ysis is presented. 
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Fig. 12. Case /j, — 2. Density and velocity fields in the xz plane. 



4.2.2. Detailed analysis 

The flow features described above tend to suggest that the 
outflows we found in this model are magneto-centrifugally 
driven. This type of outflow motion has been studied by 
many authors using self-similar techniques (Blandford & 
Payne 1982, Pelletier & Pudritz 1992) and assuming sta- 
tionarity and axisymmetry. Therefore, in order to get the 
late phase evolution of the outflow for which it is ex- 
pected that stationarity has been reached, we again use 
the model with a reduced critical density as its dynami- 
cal evolution is faster (larger timesteps in the thermally 



supported core), therefore allowing to get more easily the 
stationary regime. All figures were obtained after making 
an azimuthal average of the variables around the vertical 
axis. 

The density and velocity fields are shown in the upper 
panel of Fig. 14. The flow is similar to that obtained with a 
critical density equal to p c . To study quantitatively various 
outflows quantities, we focus on the parts of the outflow 
which are close to the equatorial plane {z/Rq < 0.32). 
Further away from the disk midplane, the outflow hits 
the inflowing material and its structure is perturbed. The 
poloidal magnetic field lines in this inner region are repre- 
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Fig. 13. Structure of the azimuthally averaged magnetic field 
in the model having fi = 2 at time t — 1.5304r//. The solid 
lines displays the poloidal magnetic field lines. They are over- 
plotted on a snapshot of the toroidal magnetic field strength. 



sented in Fig. 14 with dotted lines. The outflow properties 
are computed along one such field line, represented using 
the thick solid line in Fig. 14. One of the predictions of the 
theories mentioned above is that poloidal velocities v p and 
magnetic field B p are aligned when the outflow is in steady 
state. We plot in the upper panel of Fig. 15 the variation 
of the angle 9 they make as one moves along that selected 
field line. Apart from the very inner part of the outflow 
(z < 0.08, which corresponds to the outflow launching re- 
gion), 9 is everywhere smaller than 10 degrees, indicating 
a good alignment between the velocity and the magnetic 
field. In general, over the entire outflow region, we found 
that this angle is always smaller than 25 degrees. This is a 
good indication that the outflow has come close to reach- 
ing steady state, which is in agreement with visual inspec- 
tions of animations of this simulation. The middle panel of 
Fig. 15 gives an insight into the launching mechanism, by 
plotting the profile of the forces acting on the fluid along 
the same field line. The solid line shows the variation of 
the component of the Lorentz force along that field line. 
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Fig. 14. Upper panel: density and velocity field in the model 
having fi = 2 and a critical density p c /10 at time t = 1.67t//. 
Lower panel: structure of the azimuthally averaged poloidal 
magnetic field lines (dashed line) in the region of the outflow. 
The dashed lines show the global structure of those lines, while 
the thick line mark the selected magnetic field line along which 
some quantities will be plotted in Fig. 15. Note the different 
scale of the plot compared to Fig. 13. 



It is compared with the pressure gradient, also projected 
in the same direction. The former is clearly larger than 
the later, by one or two orders of magnitude: the outflow 
is magnetically (as opposed to thermally) driven. Finally, 
we also give the profile of the outflowing velocity along 
the magnetic field line (bottom plot of fig. 15). Because of 
the magnetic force, it increases steadily in the outflow to 
reach values of the order of 1.5 km/s. 
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Fig. 15. First panel: angle 9 (in degrees) between the poloidal 
fluid velocity and the poloidal magnetic field along the par- 
ticular poloidal field line shown in figure 14. Note that 6 is 
smaller than 10 degrees, showing a good alignment between 
both vector. Second panel: Lorentz force (solid line) and pres- 
sure force (dashed line) exerted on the fluid element along the 
same magnetic field line. At all positions, the former exceeds 
the latter by one or two orders of magnitude, indicating that 
the outflow is largely magnetically driven. Third panel: Fluid 
velocity along the magnetic field line shown in Fig. 14. The 
gas is seen to be constantly accelerated because of the action 
of the Lorentz force (see fig. 15). It reaches a maximum outflow 
velocity of the order of 1.5 km/s. 



Another important prediction of the analytical self- 
similar model (Blandford & Payne 1982) is that the angle 
between the magnetic field lines close to the disk and the 
z-axis, should be larger than 30 degrees. In Fig. 16, we 
show this angle as a function of the radius. It has been 
measured at the disk surface, defined at each radius as 
being the altitude at which the radial fluid velocity van- 
ishes. It is seen that this angle is indeed always larger than 
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Fig. 17. Ratio of ejection over accretion mass rate as a function 
of time. Upper panel is /i — 20 and lower panel /i — 2. Accretion 
and ejection mass rates are estimated on a sphere of radius R s . 
Solid lines are for R s /Ro = 0.2, dotted lines for R s /Ro = 0.4 
whereas dashed lines are for R s /Ro = 0.6. 



30 degrees except in the very center and in the outer part. 
In these two regions, no outflow is occurring as can be 
seen in Fig. 14. 
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Fig. 18. Ratio of ejection over accretion angular momentum 
rate as a function of time. Upper panel is fi = 20 and lower 
panel /i = 2. Accretion and ejection angular momentum rates 
are estimated on a sphere of radius R s . Solid lines are for 
Ra/Ro = 0.2, dotted lines for R s /Ro = 0.4 whereas dashed 
lines are for R s /Ro = 0.6. 

4.3. Mass and angular momentum fluxes 

We now present quantities which characterize globally the 
evolution of the whole accretion-ejection structure with 
time. For this purpose we again use the simulations with 
the critical density, p c /10, since they allow to follow the 
cloud evolution further. 

Figures 17 and 18 display the ratio of ejected over ac- 
creted mass and angular momentum fluxes. They are es- 
timated on spheres of various radius namely R s /Rq = 
0.2 (solid lines), R s /R = 0.4 (dotted lines) and R s /Rq = 
0.6 (dashed lines). Note that for /j, = 20, the first value of 
R s is inside the magnetic tower whereas the 2 other values 
correspond to radius higher than the equatorial radius of 
the magnetic tower. 

For /j, = 20 and R s /Ro = 0.2, the ratio of ejection 
over accretion mass rate vanishes before t = 1.3r//, then 
increases until a value of about 3-4. At this point quasi- 
stationarity is reached. This indicates that because of the 
centrifugal barrier, the gas first piles up in the inner region. 
Then the magnetic tower and the outflow described pre- 
viously extract efficiently the mass at a rate higher than 
the accretion mass rate. For the two larger values of R s , 



the ratio of ejection over accretion mass rate is smaller by 
a factor of about 3. This is due to higher accretion rates 
in the collapsing envelope than in the inner ccntrifugally 
supported structure. 

The behaviour for ji = 2, is much different. The ratio 
of ejected over accreted mass rate varies between 0.1 and 
0.6 and is therefore always smaller than 1. It increases 
with R s . In that case, the gas falls directly in the centre 
without piling up in a ccntrifugally supported structure. 
The saturated ratios obtained for the two smallest R s are 
reminiscent of typical values quoted in the literature. 

The ratio of ejection over accretion angular momen- 
tum rates have a similar behaviour than the ratio of ejec- 
tion over accretion mass rates. However, we note that for 
fi = 20, the former is smaller than the later by a factor 1 
to 2 whereas for /j, = 2, the contrary is true (the ratio be- 
ing as high as 3). These differences are due to the fact that 
in case fi = 20, the transportation of the angular momen- 
tum is weak since the magnetic tension is weak. Therefore 
the angular momentum is not transferred efficiently and 
is mostly advected with the gas. Since the gas which is 
accreted comes from larger radius than the gas which is 
ejected, the later has on average a larger angular momen- 
tum than the former. On the contrary, in case /i = 2, the 
gas is efficiently braked near the equatorial plane whereas 
it is azimuthally accelerated at higher altitude therefore 
carrying with it a higher angular momentum. 

5. Comparison with observations 

Here we qualitatively discuss comparisons between the 
models presented in the previous sections and various ob- 
servations. One of the difficulties in carrying out detailed 
comparisons between observations and models of proto- 
stellar dense core is the need for sources sufficiently con- 
strained observationally. 

5.1. The case of IRAM04191 

In this respect, the 1.5 solar mass, young class source 
IRAM04191 (Andre et al. 1999, Belloche et al. 2002) lo- 
cated in the Taurus molecular cloud, is a nice example. 
In this elongated source, an outflow perpendicular to the 
major axis of the source has been observed suggesting that 
the rotation axis is also perpendicular to the major axis. 
With these assumptions, the rotation velocity has been 
measured. Moreover, the radial velocities and the column 
density profile are known as well from radiative modeling 
of the line profiles. A dynamical age of ~2 10 4 years has 
been estimated from the characteristic of the flow. Finally 
no disk has been detected in this source, the upper limit 
for the disk radius being around 15 AU. 

Various attempts have been made to compare these 
profiles with hydrodynamical models, starting initially 
with critical Bonnor-Ebcrt sphere in rotation (Belloche 
2002, Hennebelle et al. 2004b, Lcsaffrc et al. 2005). These 
models fail to reproduce IRAM04191 for the following rea- 
sons. First, the infall velocity (~ 0.15 km/s) is too large 
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at r ~ 2000 AU in the model (0.2-0.3 km/s). Second the 
column density in the inner part (r < 1000 AU) is too 
large in the model. In the model, the large column den- 
sity in the inner region is due to the influence of rotation 
(see Fig. 1) which provides a support to the infalling gas. 
Self-consistently, the rotation curve cannot be reproduced 
under simple assumptions on the initial angular momen- 
tum distribution. In particular, the rotation velocity of 
the hydrodynamical model tends to be too high in the in- 
ner part. Another important related disagreement, is the 
absence of a big massive disk in IRAM04191 like the one 
produced in the hydrodynamical simulation. 

Although no detailed comparison has been carried out, 
it seems worthwhile to compare with the magnetized mod- 
els presented in this paper qualitatively. The first interest- 
ing aspect is the infall velocity which is smaller in the outer 
part than in the hydrodynamical case because of magnetic 
support. Although in the models presented here, it is still 
too large to reproduce the infall of IRAM04191, starting 
with near equilibrium configuration will help to reduce 
it further. However, running specific models dedicated to 
this particular case is beyond the scope of this paper. The 
second aspect is the rotation curve, which is much flatter 
in the fi = 2 — 5 cases than in the hydrodynamical case. 
This is also in better agreement with the rotation curve 
inferred for IRAM04191. Finally, the models with [i < 5 
do not show the presence of hundred AU size disk unlike 
the hydrodynamical model. It is indeed, extremely diffi- 
cult to reconcile the absence of disk and the presence of 
rotation within the framework of hydrodynamical models. 

Another interesting aspect is the outflow which is ob- 
served in IRAM04191 (Andre et al. 1999). Its velocity is 
about 5-10 km/s. This is qualitatively in good agreement 
with the outflows spontaneously launched in the MHD 
models. Quantitatively, however, this is 2 to 5 times larger 
than the values obtained in this paper (note that outflow 
velocities up to ~3 km s^ 1 are produced in our simu- 
lations). Nevertheless, as recalled previously, one should 
remember that the speed of the outflow is related to the 
rotation velocity achieved by the fluid particles. Since the 
physics of the first Larson core and the second collapse 
are not treated in this work, further contraction of the 
gas is prevented and therefore the velocity of the outflow 
is reduced. It seems therefore reasonable to assume that 
a better treatment of the first Larson core will yield to 
faster outflows (see Banerjee & Pudritz 2006, Machida et 
al. 2007). 

5.2. Observations of class sources with disk 

While disks are commonly observed during the late stage 
of star formation, i.e. class I, class II or T Taury phase (e.g. 
Watson et al. 2007), disks are more difficult to observe dur- 
ing the class phase and therefore much less constrained 
(Mundy et al. 2000). Nevertheless, various studies report 
class protostars having disks of masses between 1 and 
10 percents of their envelope masses (Looney et al. 2003, 



Jorgcnsen et al. 2007) giving a typical mass of about 0.1 
solar mass. 

Since the age of these sources is not well known and the 
density and velocity profiles are not available, it appears 
difficult to reach solid conclusions. This may nevertheless 
indicate that in these sources, the magnetic field is not too 
strong with typically fi > 5. 

6. Conclusion 

Using RAMSES, we performed 3D simulations of a mag- 
netised collapsing dense core. We explored the effect of 
the initial magnetic field strength varying the value of the 
mass-to-flux over critical mass-to-flux ratio, /x, from 1000 
to 2. The cloud evolution is significantly modified for val- 
ues as large as [i = 20. This is due to the strong am- 
plification of the radial and azimuthal components of the 
magnetic field induced by the differential motions arising 
in the collapsing cloud. 

We also developed semi-analytical models that predict 
some of the core envelope properties and compared them 
with the simulations results showing reasonable agree- 
ments. 

For [i — 20, we find that magnetic braking is negligi- 
ble and that consequently a centrifugally supported disk 
forms. A magnetic tower, generated by the twisting of the 
field lines, forms and expands reducing the mass of the 
centrifugally supported structures. A faster outflow is then 
triggered from the thermally supported central core. 

For [i smaller than 5, no centrifugally supported disk 
forms for two reasons. First the collapse occurs primarily 
along the field lines which means that less angular mo- 
mentum is delivered to the inner parts. Second, strong 
magnetic braking extracts the angular momentum from 
the disk. The question as to whether a disk will form 
at later stage remains nevertheless open. In addition, an 
outflow is triggered from the thermally supported core in 
that case. Detailed investigations have been performed in 
the case /i = 2. They reveal that the outflow reaches a 
quasi-steady state regime and features many characteris- 
tics of the magneto-centrifugally driven outflow models 
studied analytically in the literature (e.g. Blandford & 
Payne 1982, Pudritz et al. 2006). 

In a companion paper, we study the fragmentation of 
the collapsing dense core paying particular attention to 
the strength of the magnetic field. The analysis developed 
in the present paper is then used to interpret the numerical 
results obtained in the context of fragmentation. 
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Appendix A: An analytical model for the 
expansion of the magnetic tower with 
self-gravity 

We present an analytical model to investigate the mechanism 
responsible for the expansion of the self-gravitating disk. 

In the simulation, the problem appears to be axisymmetric 
making it bidimentional. It is also evidently time dependent. 
For the purpose of reducing the complexity, we consider a self- 
gravitating layer which varies along the z-axis only, instead of 
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an axisymmetric structure. We therefore replace, the azimuthal 
fields, Bg and V$ by the plane-parallel transverse fields B x (z) 
and V x {z). It is initially threaded by a vertical and uniform 
magnetic field, B®. The transverse velocity, V x (z), initially gen- 
erates a transverse magnetic field, B x (z) which modifies the 
vertical equilibrium. By doing so, we ignore the effect of the 
curvature inherent to the axisymmetric structure. 

For simplicity, we make the approximation that the layer is 
in mechanical equilibrium and use Eq. (9) written previously. 
Let <T = 2j pdz be the column density. With the Poisson equa- 
tion, we get 

d z tp = -2%Ga (A.l) 
Thus, defining 

p = p/po , a = a/V2p cI/ttG, B x = B x /B , (A.2) 
where Bo = \f^Poc 2 3 , we get, since <r(0) = and B x (0) = 0, 

~2 

p + a 2 + ^ = l, (A.3) 

where po is the density at z = 0. 

In order to compute B x , we write the transverse momentum 
and induction equations using the Lagrangian coordinates (Shu 
1983). With t — t/to where to = l/\/27rpoG, and z — z/zo 
where z — a /2p we have 

dfV x = B° z dfB x , (A A) 



and 




(A.5) 



This equation shows that, in the model, the growth of the 
toroidal field is triggered only by the vertical gradient of V x , 
the transverse velocity. It is worth stressing that this equation 
does not include all the physical relevant terms in the prob- 
lem. In particular, in the simulation the growth of the toroidal 
field is largely due to the twisting of the radial component of 
the magnetic field by the differential rotation which cannot be 
included in a plane parallel geometry. We note however that 
Fig. 11 clearly shows vertical gradients of Ve (z > 1CP 4 pc at 
time 1.116 and 1.125 r//). 

Equations (A.3), (A. 4) and (A.5) are to be solved. Despite 
the simplifications, i.e. dependence on z only and mechanical 
equilibrium in the vertical direction, they are still complex two 
variable non-linear equations. Since the disk is symmetric with 
respect to the equatorial plane, the boundary conditions are 
B x (0) = and d a V x (0) = 0. 

In order to illustrate the origin of the expansion due to the 
growth of the toroidal magnetic component, we consider as 
initial conditions a self-gravitating layer with a vanishing B x 
at time t = 0. In that case the solution is simply p = 1 — a 2 . 
In term of the z variable, this writes p(z) = l/ch(2z) 2 (Spitzer 
1942, Ledoux 1951, Curry 2000). 

Since obtaining exact solutions of Eqs. (A.3)-(A.5) seems 
to be difficult, we seek approximated solutions of the equations 
written above. To this purpose, we replace Eq. (A.5) by 

d T B x = (l-a 2 )d 3 V x , (A.6) 

that is to say, we assume that the density in Eq. (A.5), is the 
density of the unmagnetised solution. Strictly speaking, this 



approximation holds as long as the density has not been mod- 
ified by the magnetic field significantly. With this assumption, 
the time and space variables separate, making it easy to find 
solutions as 

V x (t, a) = V cos(V6B°t) x (1 - 3a 2 ), (A.7) 

B x (t, a) = -VeVo sm(V6B°t) x (5 - a 3 ). (A.8) 

Using Eq. (A.3), we obtain the density as a function of <r and 
t 

p = 1 - a 2 - SV() 2 sin 2 (V6B°t) (a - a'') 2 . (A.9) 

using the relation da — 2pdz, it is possible to obtain z as a 
function of a and p. 

Figure A.l shows p and B x as given by Eqs. (A.8) and (A.9) 
at four times for Vq = 2 and B° z = 1. The behaviour is very 
similar to the evolution displayed in Fig. 11. As B x is growing, 
the layer is expanding because of the magnetic pressure. 



